A Scale-free Network with Boolean Dynamics as a Function of Connectivity 
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I. INTRODUCTION 



This paper blends two promising concepts of complex 
systems, namely the boolean networks [1,2] and models 
of growing networks [3,4]. The former was introduced 
by Kauffman [1,2] in 1969. One important thing about 
Kauffman's model is that it was very well received not 
just by the physical but by the biological community as 
well. The simplicity and richness of behaviors of this 
model created a huge field of research, with many types 
of different boolean networks [5]. Moreover, this model 
became a good candidate to explain biological problems 
such as cell differentiation, gene expression, protein inter- 
action and genetic regulatory networks [6-8]. A boolean 
network is a complex dynamical system constituted of 
logical variables connected by logical functions. The sim- 
plest boolean model has two parameters, the number of 
logical variables N and the number of inputs of their 
boolean functions K, where K varies from K = 1 (the 
function has just one input) to K — N (all variables con- 
nected with all the others including itself). The boolean 
functions are chosen randomly in the beginning and they 
are kept during the dynamics. The network updates syn- 
chronously, meaning tha all nodes refresh their state at 
the same time. This kind of network exhibit an ordered 
dynamic for K — 1 and a chaotic one when K = N . 
However, the network self-organizes showing an interest- 
ing level of order and complexity at K = 2, the edge of 
chaos. In addition to Darwin's natural selection and ran- 
dom mutation, Kauffman's idea is that self-organization 
and random dynamic can be responsible for the complex- 
ity observed in nature. 

The models of growing networks introduced by 
Barabasi and Albert [3,4] are in the other side. In these 
works was proposed a very elegant way of creating a self- 
organized scale-free structure. In contrast to most types 
of networks previously known , the scale- free ones exhibit 
a power-law distribution of connections, P{k) ~ k^ 1 
where 7 is the so called scale- free exponent [9]. This 
behavior is similar to what happens with lengths in frac- 



tal structures. The mechanism that makes possible the 
arising of such property is the so called preferential at- 
tachment, that is, the probability of an old node to receive 
a new link from a new node is proportional to the num- 
ber of pre-existent connections in this node. The growing 
nature of these networks allows the study of a wide type 
of networks, such as the World Wide Web, neural, social 
relations, disease spreading, voting, citations of scientific 
papers, movie's actors interconnections, etc [10-21]. The 
second characteristic of a scale-free network is the small 
world effect, which means that the shortest distance be- 
tween two nodes is of the order of ln(iV), where N is the 
size of the network. In other words, no matter how large 
the network could be, any two nodes are connected via 
a relatively small chain of links. For instance in WWW 
any document can be reached with less of 19 mouse clicks. 
The third feature of scale- free graphs is the very good tol- 
erance to random removal of a significant fraction of its 
nodes conjugated with a high vulnerability to directional 
attacks to the most connected nodes (hubs). 

Although these two branches of research are in plenty 
expansion, with a huge amount of publications, the com- 
bination of both are just in the beginning, and little is 
known about the behavior of the random networks under 
scale-free topology. 

Since scale-free networks describe more realistically the 
interactions among members of any types of network, the 
boolean dynamic can help to understand how some sort 
of dynamics flow inside those groups with that topology. 
From the biologic point of view, the introduction of the 
scale- free structure can be a decisive step to describe in a 
qualitative and quantitative fashion what is observed in 
genetic networks, metabolic pathways and protein inter- 
action networks. All these networks have developed from 
very simple ones under the preassure of natural selection. 
They have evolved by gradual changes (mutations) that, 
simultaneously keep its functionality. Therefore we ex- 
pect that scale-free networks are a more realistic descrip- 
tion since they are growing graphs. It is worth mention- 
ing, that scale-free network are more robust to random 
errors than homogeneous random models. Experimental 
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evidences have been found in cellular networks of some 
living beings, as example the yeast Saccharomyces cere- 
visiae [23]. 

Few papers have been published using boolean dy- 
namics in scale-free networks. Fox and Hill used a ran- 
dom dynamics in a network with maximum connectivity 
K max = 30 to simulate the regulation of gene expression 
[22]. Aldana and Cluzel analytically demonstrated the 
existence of a phase transition for values of the scale- free 
exponent in the open interval (2, 2.5) was analytically 
demonstrated in the random dynamics [24] . 

In this work we study scale-free networks with a deter- 
ministic boolean dynamics using numerical simulation. 
We consider that the dynamics is driven only by AND 
and XOR functions. These functions are controlled by 
an external parameter q. This simple dynamics allows us 
to simulate large networks. We consider the Hamming 
distance D and the number of l's M for asymptotic times 
and different values of q. After averaging for several ini- 
tial conditions we find that these quantities vary with q 
as power laws. 

In the following section we define the scale free net- 
works and the boolean dynamics. Section III present our 
results for the numerical simulations, which are then dis- 
cussed in Section IV. A brief summary is given in Section 
V. 



II. THE NETWORKS AND THE DYNAMICS 

The computer simulation was performed in two ma- 
jor parts. The first part was the growth of the network, 
which corresponds to the static part of the simulation 
since the network is unchangeable during the dynamics. 
We grew scale-free networks with minimum connectivity 
kmin — 1 and k m i n — 2 by using the Growing Network 
with Re- direction algorithm [25]. The minimum connec- 
tivity corresponds to the smallest number of links that a 
node can have. The two kinds of network are grown as 
follows: 

• k m in = 1 - a new node is linked to only one old 
node; we select an old node with uniform probabil- 
ity; then the link with the old node is established 
with probability 1 — r or it is redirected to the an- 
cestor of the old node with probability r; 

• k m in = 2 - a new node is linked to two old nodes; 
we select an old node with uniform probability; 
then one link with the old node is established with 
probability 1 — r or it is redirected to one of the 
two ancestors of that node with probability r; we 
repeat the same procedure for the other link. 

In both cases, the initial condition consists of three nodes 
with cyclic connections. This algorithm creates a scale- 
free network with 7 = + 1. For instance, when 
r = 0.5, the network has 7 = 3, corresponding to a 



growth with linear preferential attachment. We use dif- 
ferent values of r for the two kinds of networks. 

The second part of the simulation is the dynamic itself. 
Once the network is made of TV nodes connected by links, 
to each node i is assigned a logical variable ai(t). The 
state of the network at time t is represented by a set of 
boolean variables (<7i(t), 02(f), os(t), <7jv(i)). At each 
time step, the state of a node i is defined in the following 
way: 

i) for kmin = 2, ai(t + 1) is given by 
F i (ai 1 (t),ai 2 (t), ...,ai k .(t)). In other words, the state of 
a node i at t + 1 is a function of the states at t of all nodes 
linked with i, namely cr^ (t), o~i 2 (i), o~i k (t), where ki is 
the connectivity of i-th node. 

ii) for k min = 1, the state of a node is defined in a dif- 
ferent way, since there are nodes with just one link. Then 
cri(t + 1) is given by Fj(cr; (t), a ir (t), a i2 (t), a ik . (t)). 
Now, the state of a node i at t + 1 is a function of the 
state at t of all nodes linked to it and of its own state as 
well. In this way we have always a function with at least 
two inputs. 

Finally we define the function Fi(a): 

p = f AND iih<q m 
1 XXOR otherwise 1 ' 

where q is a threshold parameter that controls how the 
logical functions AND and XOR spread into the network. 
The AND and XOR logical functions were introduced 
in order to simplify the model, avoiding the necessity 
of defining 2 2 different boolean relations for each node. 
It is know from previous works that the AND function 
leads to an ordered regime (with two fixed points, where 
all variables are 0's or l's) while the XOR introduces a 
more chaotic component to the dynamic. These kinds of 
functions are necessary to study biologic networks since 
they are the boolean counterparts of real reactions in cell 
regulatory system [2,7]. Once in the scale-free network, 
we have many nodes poorly connected (small k) and few 
nodes highly connected (large k), equation (1) we can set 
a balance between chaotic and ordered dynamics inside 
the network. 

An initial state S = (<7i(0), <r 2 (0), 03(0), a N (0)) is 
created by assigning randomly 0's and l's to all nodes. 
A copy £ = (oT(0),CT 2 "(0),ai'(0), 077(0)) of the initial 
state is also created and we introduce a damage by chang- 
ing the value of one randomly chosen node. Both the ini- 
tial state and the damaged state evolve under the control 
of equation (1). Once the new state of all nodes is calcu- 
lated the entire network is updated (synchronous update) 
and the system goes to the next Monte Carlo time step 
(mcs). Note that, except for the random choice of the 
initial state, the dynamics is deterministic. This means 
that the boolean functions act with a probability p = 1 , 
and this is another simplification to the dynamics. Ran- 
dom variables are important in order to simulate real 
genetic networks under influence of many uncertainties 
like biologic variability, experimental noise and interact- 
ing variables impossible to quantify. However, this sim- 
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plification allows us to consider large networks without 
loss of complexity in the dynamical behavior. 



III. THE SIMULATION DATA AND RESULTS 



The dynamical behavior is characterized by two quan- 
tities. The first is the average density of l's, namely 



1 N 

M(q,t)= lim $>(*)> • 
N—>oo iv '—^ 
i=l 



(2) 



The second is the average of the Hamming distance, 
that is defined by 
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FIG. 1. Plots of M(q, t) vs t and D(q, t) vs t for a network 
with N = 8 x 10 4 , k min = 1 and r = 0.5. 



1 N 

D{q,t) = lim (— V 

i=l 



<Ji{t) ~ <Ji{t)\) 



(3) 



Here, (. . .} is the average for different initial conditions. 
An initial condition is given not only by the initial states 
(S, S) but also by one set of links of a grown network 
with a specific 7. 

After a very short transient (less than 30 time steps) 
these quantities reach the stationary values M(q) and 
D(q), as it is shown in Fig. 1 for a network with N = 
8 x 10^ nodes, k m i n = 1 and r = 0.5. The stationary 
values, defined as 



rt+T 

M(q) = lim l/T / M(q,t')dt' 



D{q) = lim l/T 

1 —*oo 



rt+T 



D(q,t')dt' , 



were determined by discarding the first 50 time steps and 
by making a time average until t — 200 time steps. We 
can also see in the figure that the stationary values D(q) 
and M{q) depend on q. We are interested in this depen- 
dence. It is worth mentioning that a similar behavior is 
found for networks with different values of r and also for 
all networks with k m i n = 2. 



In order to consider the finite size effects, we grew net- 
works with N = 1 x 10 4 , N = 2 x 10 4 , N = 4 x 10 4 and 



N = 8 x 10 for networks with k„ 



1 and k r . 



2. 



The sample averages were performed at least with 10 2 
samples (small q and large TV). For large q we have used 
up to 3 x 10 4 random initial conditions. 
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FIG. 2. Log-log plots of M(q) vs q and D(q) vs q for 
network with fc min = 1, r = 0.5 and different N. 
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We can see from Fig. 2 that D(q) and M(q) have the 
following asymptotic power-law behaviors: 



much larger than the ones we consider here. 
10° 



M(q) ~ q- m 
D(q) ~ q- d 



(4) 
(5) 



Moreover, we can also see in the same figure, finite 
size effects for large q by comparing the behavior of the 
smallest network (N = 1 x 10 4 ) with the largest one ( 
N = 8 x 10 4 ). In order to evaluate the exponents, we 
eliminate the points affected by finite size effect and co- 
alesce all different sets. Finally we do a best fit. This 
is shown in Fig. 3, for the case defined in Fig. 2. There, 
the first point (q = 2) was not considered in the best fit. 
Note that we have evaluated the exponents by consider- 
ing approximately 2 orders of magnitude in the q variable 
and that the fit is very good. In fact, in all fitted data, 
we obtained a correlation coefficient larger that 0.999. 
10 




FIG. 3. Log-log plots of M(q) vs q and D(q) vs q for a net- 
work defined in Fig. 2. It shows the best fit of the coalesced 
sets. 

Plots of D(q) versus q are shown in Fig. 4 for networks 
with kmin = 1 and different values of r ((a)r = 0.35 and 
(b)r = 0.80). Note that the finite size effects are present 
for q w 260 for the r = 0.8 and N = 8 x 10 4 case. This 
implies that finite size effects appear for q « k* , where 
k* is the static cut-off. In fact, when we are evaluating 
the connectivity distribution P(k) for this case, we expect 
that finite size effects appear when k ~ k* = N r « 8x 10 3 
[26]. In a numerical simulation of the exponent 7 for 
this case, we found that the finite size effects appear for 
k ~ 1000. Therefore, if we want to evaluate the expo- 
nents with 3 orders of magnitude, we must grow networks 
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FIG. 4. Log-log plots of D(q) vs q for a network with 
kmin = 1 and (a)r = 0.35 and (b)r = 0.80. The best fits 
of the coalesced sets, after the elimination of the points af- 
fected by finite size and the discard of the first point, furnish 
(a)d = 3.98 and (b)d = 1.59. 

We follow this procedure to calculate all the other ex- 
ponents. They are displayed in table I. The table also 
show the values of the scaling-free exponent evaluated 
numerically j n and the exact value 7 = 1 + r _1 . In the 
numerically evaluation of 7„ we considered the same net- 
works used in the dynamics with 10 5 samples. We can 
see that the exponents m and d change with r and k m i n . 
It is worth mentioning that the asymptotic value of D(q) 
is independent of the initial amount of damage. Only the 
short time behavior depends on it. In particular if the 
initial damage is larger than D(q), D(q, t) shows a decay 
to the stationary value. 

TABLE I. Exponents m and d for kmin = 1 and kmin = 2 for 
different values of r. It is also shown the scaling-free exponent 
evaluated numerically j n and the exact one 7. 





kmin — 1 


kmin — 2 






r 


m 


d 


m 


d 






0.35 


2.25 (3) 


3.98 (5) 


2.77 (1) 


4.10 (2) 


3.50 


3.86 


0.5 


1.71 (1) 


2.89 (2) 


2.00 (1) 


2.80 (3) 


2.92 


3.00 


0.65 


1.38 (2) 


2.15 (5) 






2.52 


2.54 


0.8 


1.16 (3) 


1.59 (1) 


0.53 (1) 


0.70 (1) 


2.27 


2.25 
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In Fig. 5 we show the temporal evolution of D(q, t)/q 



for the network with k„ 



2 and r = 0.5. Note that 



after an initial transient, all the curves collapse. 
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FIG. 5. Temporal evolution of D(q, t) for different q and 
N, where k mirL = 2 and r = 0.5. (see legend). The y axis was 
rescaled by D(q, t)/q~ d in order to show the collapse. 



IV. DISCUSSION OF RESULTS 

In order to discuss the simulation results, we present a 
simple approximation for the evaluation of the exponents 
m and d for the case k min = 1 , in which we neglect all the 
correlations between nodes. Note that argument is also 
valid for k min = 2. The dynamics defines two sublattices 
of nodes depending on the connectivity k of each node. 
The sublattices a and b consist of all nodes with k < q 
and k > q, respectively. Then the nodes of subllatice a(b) 
evolve by AND (XOR) operations. We are interested in 
the steady state, where the distribution of connectivity 
is described by P(k) ~ fc -7 . In this regime, we consider 
that a typical node i of the sublattice a with connectivity 
k has Ui = 0, because the AND operation of the nodes 
o-j, <7i . . . (Tfc is if at least one of these k + 1 nodes is 
zero. On the other hand, a node of the sublattice b has 
a probability 1/2 of being nonzero, because half of all 
possible configurations of ai, a\ . . . ah, under XOR oper- 
ation, give Gi — 1. The average density can be written in 
terms of the sublattices as M(q) = M a {q) + M b {q) ■ Since 
M a (q) w 0, we have 



M(q) « M b 



1 f c 



P(k) dk ~ q 



'(7-1) 



(6) 



implying that m = 7 — 1. To evaluate D(q) we must 
consider a lattice and its copy that have evolved from 
two different initial states. Again we have that D{q) = 
D a (q) + Dt(q), where Dj(q) is the contribution of the 
sublattice j(j = a,b). So, we need the number of nodes 
that have different values in the lattice and its copy for 
each one of the sublattices. Since the nodes of the sublat- 
tice a have the value in both lattices, D a (q) rts 0. Only 
the nodes of the sublattice b can be different from 0. The 
probability that the same node be 1 at one lattice and 
in its copy is 1/4, if we neglect all correlations. Therefore 
we can write that 



D(q) « D b (q) « \ J°° P(k) dk ~ q 



(7-1) 



(7) 



Using the definition of the exponent d, we find that 
d = m = 7 — 1. From Table I, we see that the values 
of the exponents are different from the ones predicted 
by our approximated evaluation. This suggests that cor- 
relations are important in this problem. To obtain a 
better understanding we studied the steady-state behav- 
ior of the sublattices for networks with N = 10 4 nodes. 
The results of the sublattice a are shown in Table II in 
percentages of the quantities characterizing the complete 
network. 

Let us first analyze the case with fe mi „ = 1. From Ta- 
ble II, we note that the sublattice a is irrelevant for the 
Hamming distance independently of the value of r. 

TABLE II. Sublattice percents of M a (q) and D a (q) for 
q — 5, 10 and 16 with k mi „ = 1, N = 10 4 and different values 
of r. 





M a 


Da \ 


q 1 r 


0.35 


0.5 


0.8 


0.35 


0.5 


0.8 


5 


7% 


6% 


4% 


0.03% 


0.03% 


0.01% 


10 


18% 


16% 


9% 


1% 


0.3% 


0.01% 


16 


25% 


24% 


15% 


4% 


1.3% 


0.3% 
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On the other hand, the contribution of this sublattice 
for the average density is relevant. This contribution for 
r = 0.35 is larger than the one for r = 0.8. This indi- 
cates that the sublattice a is almost irrelevant for r = 1. 
To get a deep insight we analized M b {q, k) and D b (q, k), 
respectively, the average density and the average Ham- 
ming distance of nodes with connectivity k for a fixed 
parameter q. Note that both definitions involve only the 
nodes of sublattice b. From the numerical simulations 
we obtain that the behavior of these quantities can be 
described by 

M b (q,k) ~A q k- a , (8) 

D b (q,k) ~ B g k~P , (9) 

with a and (3 having values approximately equal to 7„, 
the numerical estimate of 7. This implies that a = f3 = 7. 
Therefore, the sublattice quantities can be expressed as 

M b (q) = f™M(q,k)dk~A q q-( a -V , (10) 

D b {q) = f~ D{q,k) dk~ B q q-W-V . (11) 

In Fig. 6 it is shown the graph M b (q, k)xk for k min = 1, 
r = 0.5, TV = 10 4 and q = 5, 10, 16. A best fit furnishes 
a = 2.89 = 7„ and an intercept A q , both independent 
of q. Similar results are obtained for the other values of 
the parameter r. If we neglect M a , then we can write 
that M(q) ~ M b and m = a- l = j-l. Then, the 
exponent m should have a value near 7 — 1. This can 
be verified in Table I and the difference between 7 — 1 
and to, which is around 20% for r = 0.35 and 7% for 
r = 0.8, is due mainly to the contribution of the sublat- 
tice a. Since this contribution becomes smaller and the 
exponent to approaches 7 — 1 as r — ► 1, we conjecture 
that to = 7 — 1 = 1 when r — 1. 
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FIG. 6. Log-lo plots of the sublattice b quantities 
Mb{q, k) vs k and Db(q, k) vs k for different q, where k min — 1, 
r = 0.5 and N = 10 4 . 

Fig. 6 displays D b (q,k) x k. A best fit furnishes 
(3 = 2.89 w 7„ independent of q and an intercept B q that 
depends on q. This dependence is due to the correlation 
between the nodes of the two sublattices. Similar results 
are obtained for different values of r. Since D a w in- 
dependently of r, we have that D{q) w D b (q) <~ Bqq 13 ^ 1 , 
with (3—1 = 7 — 1. Then the exponent d should have 
a value different from 7 — 1 because of B q . This is con- 
firmed in the simulations, with the intercept presenting 
a power law behavior B q ~ q^^ 1 . For r = 0.5, (3\ w 1 
and we have that d w 7. When r becomes larger, the 
straight lines for different g in a log-log plot become very 
close, implying that (3\ — > 0. Moreover, the exponent d 
approaches 7 — 1 as r — > 1. Then, we conjecture that 
d — 7 — 1 = 1 for r = 1 . 

It is worth mentioning that for r = 1, a new node is 
always redirected to the ancestor node. It means that we 
have three hubs, the initial nodes, with all other nodes 
connected to them. In fact we have a largest hub, with 
61% of the all nodes connected to them, a smallest one 
(11% of connections) and a third hub with 28% of all 
links. During the dynamics, the XOR operation is ap- 
plied to the hubs and AND to the nodes with a single 
link. 
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FIG. 7. Exponents d and m — 0.2 as a function of r for 

kmin = 1- 

The exponents d and m — 0.2 as a function of r for 
kmin = 1 are shown in Fig. 7. Note that we included the 
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values m = d = 1 for r = 1 in the data. Although we 
have only 5 points, we can do a best fit. We obtain that 
d = 0.963 - 2.836 ln(r) and m - 0.2 = 0.795r-°- 908 with 
correlation coefficients around —0.9995. 

Now, let us discuss the case k m i n = 2. We can see 
from Table 111 that the sublattice a is now relevant. In 
fact, as long as r — > 1, this sublattice is more important 
than subllattice b. Moreover, the nodes with connectivity 
k = 2 are responsible for the major contribution of D(q) 
or M(q). For example, when we have that r = 0.8 and 
q = 10, they correspond to 88% of D and to 83% of M. 
The scenario for r = 1 is as follows: sublattice b has no 
direct contribution (D^ w and M& « 0) and only the 
nodes with k = k m i n = 2 are responsible for D(q) ^ 
and M(q) ^ 0. Now, we have three hubs, the initial 
nodes, with all other nodes connected to them by two 
links. The largest hub and the smallest one have, respec- 
tively, 52% and 17% of the all nodes connected to them. 
During the dynamics, the XOR operation is applied to 
the hubs and AND to the other nodes. These arguments 
imply that the exponents d and m are not related to 7 
in the same way as the previous case (k min = 1). 



TABLE III. Sublattice percents of M a (q) and D a (q) for 
q = 10 and 20 or 25 with k min = 2, N = 10 4 and different 
values of r. 





M a 




q 1 r 


0.35 


0.5 


0.8 


0.35 


0.5 


0.8 


10 


24% 


43% 


84% 


32% 


53% 


89% 


20 


17% 


43% 




22% 


53% 




25 






92% 






95% 
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V. SUMMARY 

We have studied a deterministic boolean dynamics 
with two boolean functions (AND and XOR) controlled 
by an external parameter q. We have considered two 
distinct networks with minimum connectivity given by 
kmin = 1 and k m i n = 2. In the first case, the state of a 
node at time t + 1 is a function of all connected nodes 
plus its own value at previous time t. In the second case, 
the state of a node at t + 1 depends only on the states 
of all connected nodes at time t. We have grown net- 
works with different scale-free exponents by changing a 
parameter r. The finite size effect were take into account 
by considering networks with different sizes. We have 
shown that the density of the Hamming distance D and 
the density of l's M as a function of q have a power law 
tail for asymptotic times. 

It seems that the exponents d and m characterizing the 
behaviors of D and M, respectively, depend on k m i n and 
r, this means that the exponents depend on the details 
of the dynamics. 
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